% Trade Elasticity: computed as response to a 100bps decline in US tariff in BGP

clear all; clear global; clc; close all;

load('../Set_postestimation_use.mat')

tau_A = 0.053;
tau_B = 0.038;

tau_bA = 0.053;
tau_bB = 0.041;

Par.tau_A = tau_A;
Par.tau_B = tau_B;

Par.L_A   = 1;
Par.L_B   = 1;

Par.time_ss  = 2*Par.time_ss;
Par.time_81  = Par.time1;
Par.time_95  = 20/Par.dt;
Par.time_add = (200+Par.time_81*Par.dt)/Par.dt-Par.time_ss+1;

Par.num_r = Par.num_r/2;


horizons = [1:5 10:10:150]; 
trfA_grid = Par.trf_A-0.01;

[Par,Step,Trans] = Trade_ctrf(Par);

rev_b   = 1; % Revenue back to consumer
retal_b = 0; % 1 if foreign retaliation

Par.rev_US = rev_b;
Par.rev_FN = rev_b;

% ============================================
%                 SOLUTIONS
% ============================================

Param = structvars(Par);
for ctr0 = 1:size(Param,1); eval(Param(ctr0,:)); end

steps      = Step{1};
Trans_innov = Step{2};
Trans_fail = Step{3};
Trans_exo  = Step{4};

%% Calibrated Economy Path in BGP
% =================================================================

Tol_rr = 1*10^-6;
save Tol_rr Tol_rr

%---------------- Path in BGP

mu_init = y_init;
Q_init  = y_init;

Par.if_long = 0;

[rr_sqn, X_sqn, V_sqn, Xe_sqn, Q_sqn, Ql_sqn, mu_share, ind_neg, ...
                            prft_fnc_sqn, wage_dom_sqn, mcuts] = ...
                                Path_ctrf(mu_init,Q_init,Par,Step,Trans); % initial run

mu_init = mu_share(:,13*10^3);
Q_init  = Q_sqn(1:size(Q_sqn,1)/2,13*10^3);

[rr_sqn, X_sqn, V_sqn, Xe_sqn, Q_sqn, Ql_sqn, mu_share, ind_neg, ...
                            prft_fnc_sqn, wage_dom_sqn, mcuts] = ...
                                Path_ctrf(mu_init,Q_init,Par,Step,Trans); % extended in BGP
                            
%---------------- Variables

% Profits
% --------------------------------------------------------

prft_fnc_A = prft_fnc_sqn(1:size(prft_fnc_sqn,1)/2,:);
prft_fnc_B = prft_fnc_sqn(size(prft_fnc_sqn,1)/2+1:end,:);

wage_dom_A = wage_dom_sqn(1:size(wage_dom_sqn,1)/2,:);
wage_dom_B = wage_dom_sqn(size(wage_dom_sqn,1)/2+1:end,:);

% Qualities
% --------------------------------------------------------

QA = Q_sqn(1:size(Q_sqn,1)/2,:);
QB = Q_sqn(size(Q_sqn,1)/2+1:end,:);

QA_long = Ql_sqn(1:size(Q_sqn,1)/2,:);
QB_long = Ql_sqn(size(Q_sqn,1)/2+1:end,:);

QwA = sum(wage_dom_A.*QA+flip(1-wage_dom_B).*QA*(1+kapa)/((1+kapa)*(1+trf_B))^(1/bet));
QwB = sum(wage_dom_B.*QB+flip(1-wage_dom_A).*QB*(1+kapa)/((1+kapa)*(1+trf_A))^(1/bet));

qq_sqn = [sum(QA)./QwA; sum(QB)./QwB];
                
wq_sqn = [(1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(1,:).^(-Par.bet);
            (1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(2,:).^(-Par.bet)];  

xk_sqn = [(Par.pi_frn/Par.bet)^(1/(1-Par.bet))./wq_sqn(1,:).^(1/Par.bet)*Par.L_B;
         (Par.pi_frn_toUS/Par.bet)^(1/(1-Par.bet))./wq_sqn(2,:).^(1/Par.bet)*Par.L_A];

% Shares
% --------------------------------------------------------

MU      = mu_share(:,1:time_ss+1);
MU_long = mu_share;

McutX  = mcuts(1,1:time_ss);
McutM  = mcuts(2,1:time_ss);

% Trade
% --------------------------------------------------------

X_A = pi_frn*L_B*sum(flip(1-wage_dom_B).*QA)/bet.*wq_sqn(1,:).^(1-1/bet);
M_A = pi_frn_toUS*L_A*sum(flip(1-wage_dom_A).*QB)/bet.*wq_sqn(2,:).^(1-1/bet);
D_A = pi_dom*L_A*sum(wage_dom_A.*QA)/bet.*wq_sqn(1,:).^(1-1/bet);

IMD_A = (1+Par.trf_A)*M_A./D_A;
  

%% Structural Change: Taxes in Both Countries
% =================================================================
%

Tol_rr = 5*10^-5;
save Tol_rr Tol_rr
          

% ------------- Initialization ------------- 

mub_init = mu_share(:,time_81+1);
Qb_init  = QA(:,time_81+1);

mu_pre = mu_share(:,1:time_81);
rr_pre = rr_sqn(:,1:time_81);
QA_pre = QA(:,1:time_81);
QB_pre = QB(:,1:time_81);
QwA_pre = QwA(:,1:time_81);
QwB_pre = QwB(:,1:time_81);
McutM_pre = McutM(1,1:time_81);
McutX_pre = McutX(1,1:time_81);

QA_post = QA(:,time_81+1:end);
QB_post = QB(:,time_81+1:end);

MCUTX = [McutX; zeros(length(trfA_grid),size(McutX,2))];
MCUTM = [McutM; zeros(length(trfA_grid),size(McutM,2))];

% parpool(4)

for ctr_k = 1:length(trfA_grid)
    
    % =============== Structural Change ===================

    Par_b = Par;
    Par_b.num_r = Par.num_r/2;

    Par_b.tau_A = tau_bA;
    Par_b.tau_B = tau_bB;

    Par_b.tau_eA = 0;
    Par_b.tau_eB = 0;

    Par_b.trf_A = trfA_grid(ctr_k); 
    Par_b.trf_B = Par.trf_B+retal_b*(trfA_grid(ctr_k)-Par.trf_A);
    
    Par_b.rev_US = rev_b;
    Par_b.rev_FN = rev_b;
    
    [Par_b,Step_b,Trans_b] = Trade_ctrf(Par_b);

    % ============================================
    %                 SOLUTIONS
    % ============================================

    pi_domb = Par_b.pi_dom;
    pi_frnb = Par_b.pi_frn;
    pi_frnb_toUS = Par_b.pi_frn_toUS;

    %---------------- Path

    [rrb_sqn, Xb_sqn, Vb_sqn, Xeb_sqn, Qb_sqn, Qlb_sqn, mub_share, ind_neg_b, ...
                                prft_fncb_sqn, wage_domb_sqn, mcutsb] = ...
                                Path_ctrf(mub_init,Qb_init,Par_b,Step_b,Trans_b);

    rr_b = [rr_pre rrb_sqn(:,1:end-time_81)];

    %---------------- Variables

    % Fundamentals
    % --------------------------------------------------------

    prft_fncb_A = prft_fncb_sqn(1:size(prft_fncb_sqn,1)/2,:);
    prft_fncb_B = prft_fncb_sqn(size(prft_fncb_sqn,1)/2+1:end,:);

    wage_domb_A = wage_domb_sqn(1:size(wage_domb_sqn,1)/2,:);
    wage_domb_B = wage_domb_sqn(size(wage_domb_sqn,1)/2+1:end,:);

    % Qualities
    % --------------------------------------------------------

    QA_b = Qb_sqn(1:size(Q_sqn,1)/2,:);
    QB_b = Qb_sqn(size(Q_sqn,1)/2+1:end,:);

    QAb_long = Qlb_sqn(1:size(Q_sqn,1)/2,:);
    QBb_long = Qlb_sqn(size(Q_sqn,1)/2+1:end,:);

    QwA_b = sum(wage_domb_A.*QA_b+flip(1-wage_domb_B).*QA_b*(1+kapa)/((1+kapa)*(1+Par_b.trf_B))^(1/bet));
    QwB_b = sum(wage_domb_B.*QB_b+flip(1-wage_domb_A).*QB_b*(1+kapa)/((1+kapa)*(1+Par_b.trf_A))^(1/bet));

    qqb_sqn = [sum(QA_b)./QwA_b; sum(QB_b)./QwB_b];
                
    wqb_sqn = [(1-Par_b.bet)*Par_b.eta^(Par_b.bet-1)*qqb_sqn(1,:).^(-Par_b.bet);
                (1-Par_b.bet)*Par_b.eta^(Par_b.bet-1)*qqb_sqn(2,:).^(-Par_b.bet)];  

    xkb_sqn = [(Par_b.pi_frn/Par_b.bet)^(1/(1-Par_b.bet))./wqb_sqn(1,:).^(1/Par_b.bet)*Par.L_B;
                (Par_b.pi_frn_toUS/Par_b.bet)^(1/(1-Par_b.bet))./wqb_sqn(2,:).^(1/Par_b.bet)*Par.L_A];

    QA_orgb = [QA_pre QA_b(:,1:end-time_81)];
    QB_orgb = [QB_pre QB_b(:,1:end-time_81)];

    QwA_orgb = [QwA_pre QwA_b(:,1:end-time_81)];
    QwB_orgb = [QwB_pre QwB_b(:,1:end-time_81)];
    
    qq_orgb  = [sum(QA_orgb)./QwA_orgb; sum(QB_orgb)./QwB_orgb];

    % Shares
    % --------------------------------------------------------

    MU_b = [mu_pre mub_share(:,1:end-time_81)];

    MCUTX(ctr_k+1,:)  = [McutX_pre mcutsb(1,1:time_ss-time_81)];
    MCUTM(ctr_k+1,:)  = [McutM_pre mcutsb(2,1:time_ss-time_81)];

    % Trade elasticity
    % --------------------------------------------------------    

    X_bbA = pi_frnb*L_B*sum(flip(1-wage_domb_B).*QA_b)/bet.*wqb_sqn(1,:).^(1-1/bet);
    M_bbA = pi_frnb_toUS*L_A*sum(flip(1-wage_domb_A).*QB_b)/bet.*wqb_sqn(2,:).^(1-1/bet);
    D_bbA = pi_domb*L_A*sum(wage_domb_A.*QA_b)/bet.*wqb_sqn(1,:).^(1-1/bet);

    IMD_bbA = (1+Par_b.trf_A)*M_bbA./D_bbA;
    
    IMD_bA  = [IMD_A(:,1:time_81) IMD_bbA(:,1:end-time_81)];
        
    tr_eps_dt  = -log(IMD_bA(time_81+1)/IMD_A(time_81))/log((1+Par_b.trf_A)/(1+Par.trf_A));
    tr_eps_ann = -log(IMD_bA(time_81+1/dt)./IMD_A(time_81))/log((1+Par_b.trf_A)/(1+Par.trf_A));
    tr_eps_ss  = -log(IMD_bA(end)/IMD_A(time_81))/log((1+Par_b.trf_A)/(1+Par.trf_A));

    tr_eps_ann_long = -log(IMD_bA(time_81+1+(1:450)/dt)./IMD_A(time_81))/log((1+Par_b.trf_A)/(1+Par.trf_A));
    
end

delete(gcp('nocreate'))